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We construct a lattice kinetic scheme to study electronic flow in graphene. For this purpose, we 
first derive a basis of orthogonal polynomials, using as weight function the ultrarelativistic Fermi- 
Dirac distribution at rest. Later, we use these polynomials to expand the respective distribution in 
a moving frame, for both cases, undoped and doped graphene. In order to discretize the Boltzmann 
equation and make feasible the numerical implementation, we reduce the number of discrete points in 
momentum space to 18 by applying a Gaussian quadrature, finding that the family of representative 
wave (2+l)-vectors, that satisfies the quadrature, reconstructs a honeycomb lattice. The procedure 
and discrete model are validated by solving the Riemann problem, finding excellent agreement 
with other numerical models. In addition, we have extended the Riemann problem to the case of 
different dopings, finding that by increasing the chemical potential, the electronic fluid behaves as 
if it increases its effective viscosity. 



I. INTRODUCTION 

Since its discovery p] [2] , graphene has shown a series 
of wonderful electrical and mechanical properties, such 
as ultra-high electrical conductivity, ultra-low viscosity, 
as well as exceptional structural strength, combined with 
mechanical flexibility and optical transparence. Due to 
the special symmetries of the honeycomb lattice, elec- 
trons in graphene are shown to behave like massless 
chiral ultrarelativistic quasiparticles, propagating at a 
Fermi speed of about Vf ~ 10 6 m/s [3j [4] . This places 
graphene as an appropriate laboratory for experiments 
involving relativistic massless particles confined to a two- 
dimensional space P3- 

Electronic gas in graphene can be approached from 
a hydrodynamic perspective [BH5], behaving as a nearly 
perfect fluid reaching viscosities significantly smaller than 
those of superfluid Helium at the lambda-point. This 
has suggested the possibility of observing pre-turbulcnt 
regimes, as explicitly pointed out in Ref. [TU] and later 
confirmed by numerical simulations [TT]. All these char- 
acteristics in graphene open up the possibility of study- 
ing several phenomena known from classical fluid dy- 
namics, e.g. transport through disordered media |12j . 
Kclvin-Helmholtz and Rayleigh Benard instabilities, just 
to name a few. However, the study of these phenomena 
needs appropriate numerical tools, which take into ac- 
count both, the relativistic effects and the Fermi-Dirac 
statistics. 

Recently, a solver for relativistic fluid dynamics based 
on a minimal form of the relativistic Boltzmann equa- 
tion, whose dynamics takes place in a fully discrete 
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phase-space lattice and time, known as relativistic lat- 
tice Boltzmann (RLB), has been proposed by Mendoza et 
al. [T31 [H] (and subsequently revised in Ref. [TS] enhanc- 
ing numerical stability). This model reproduces correctly 
shock waves in quark-gluon plasmas, showing excellent 
agreement with the solution of the full Boltzmann equa- 
tion obtained by Bouras et al. using BAMPS (Boltzmann 
Approach Multi-Parton Scattering) [THl HZ] . In order to 
set up a theoretical background for the lattice version of 
the relativistic Boltzmann equation for the Boltzmann 
statistics, Romatschke et al. [TH] developed a scheme 
for an ultrarelativistic gas based on the expansion in or- 
thogonal polynomials of the Maxwell- Juttner distribu- 
tion [19] and, by following a Gauss-type quadrature pro- 
cedure, the discrete version of the distribution and the 
weight functions was calculated. This procedure was sim- 
ilar to the one used for the non-relativistic lattice Boltz- 
mann model [30H23] • This relativistic model showed very 
good agreement with theoretical data, although it was 
not compatible with a lattice, thereby requiring linear 
interpolation in the free-streaming step. Another model 
based on a quadrature procedure was developed recently 
in order to make the relativistic lattice Boltzmann model 
compatible with a lattice [23] . However, all these models 
are based on the the Maxwell- Juttner distribution, which 
is based on the Boltzmann statistics, and therefore, their 
applications to quantum systems is limited. 

In this work, we construct a family of orthogonal poly- 
nomials by using the Gram-Schmidt procedure using as 
weight function the ultrarelativistic Fermi-Dirac distri- 
bution at rest. By applying a Gauss-type quadrature, we 
find that the family of discrete (2+l)-momentum vectors, 
needed to recover the first three moments of the equilib- 
rium distribution, are fully compatible with a hexagonal 
lattice, avoiding any type of linear interpolation. This 
result is very convenient, since the crystal of graphene 
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shares the same geometry, facilitating the implementa- 
tion of boundary conditions, allowing for instance hav- 
ing a good approximation for the electronic transport 
in nanoribbons with armchair or zigzag edges [251 126] 
by implementing the typical bounce-back rule for lattice 
Boltzmann models. 

The paper is organized as follows: in Sec. [TTJ we de- 
scribe in details the expansion of the Fermi-Dirac distri- 
bution in an orthogonal basis of polynomials, and per- 
form the Gauss-type quadrature. In this section, we also 
explain the discretization procedure. In Sec. |III[ we im- 
plement the validation of our model by simulating the 
Riemann problem; and in Sec. |IV| we perform additional 
simulations for doped graphene. Finally, in Sec. |Vj we 
discuss the results and future work. 



II. MODEL DESCRIPTION 

The electronic gas in graphene can be considered as a 
gas of massless Dirac quasi-particles obeying the Fermi- 
Dirac statistics in a two-dimensional space. Thus, we de- 
fine the single-particle distribution function f(x l± ,p^) in 
phase space, being x M = (x°, x 1 , x 2 ) andp M = (p ,^ 1 ,^ 2 ) 
the time-position and energy-momentum coordinates, re- 
spectively. Here x° denotes time, x — (x l ,x 2 ) spatial co- 
ordinates, p° the energy, and p — (p 1 ,^ 2 ) the momentum 
of the particles. In the ultrarelativistic regime, we get 
p^Pfi = (in this paper we use the Einstein notation, i.e. 
repeated indexes denote summing over such indexes) . In 
our approach, we assume that the distribution function 
/ evolves according to the relativistic Boltzmann-BGK 
equation [19] , 

P^f = - f eq ) , (i) 

where r is the relaxation time, and f eq the equilibrium 
distribution, which in our case, is the relativistic Fermi- 
Dirac distribution defined by 

leqix^.P^) = e {p aU c-^/ kBT + 1 > ( 2 ) 

with T the temperature, ks the Boltzmann constant, 
the macroscopic (2+1) -velocity of the fluid [T9l 127] . 
and /j, the chemical potential. The relation between 
the Lorentz-invariant If 1 and the classical velocity u — 
(u , u 2 ) is given by C/ M = ^{vf, u , u 2 ), with vf being the 
Fermi speed and 7=1/ \Jl — u 2 /v F . 

A. Moment expansion 

Here, we perform an expansion of the Fermi-Dirac dis- 
tribution, Eq. ([2]), in an orthogonal basis of polynomi- 
als. In our case, since we are interested in the hydro- 
dynamic regime, we will truncate the expansion preserv- 
ing only the polynomials up to second order, although 



achieving higher orders is also possible by using the 
same procedure. In particular, we need to reproduce 
the first three moments of the equilibrium Fermi-Dirac 
distribution, namely (l)( eg) , (p a )( eq ), and (p a p p )( eq ) for 
a,/3 = 0,1,2. The angular brackets denote expectation 
values using the distribution / via (Q) = J d/i Qf, with 
d/i = d 2 p/2p°(2Tr) 2 , and the subscript ( e? ) indicates that 
the equilibrium distribution f eq is taken instead of /. 

This method was originally introduced by Grad [28 
who expanded the Maxwell-Boltzmann distribution in 
Hermite polynomials, based on the fact that they are 
orthogonal, using as weight function the Maxwellian dis- 
tribution at rest. In this spirit, we will derive a new 
basis of polynomials that are orthogonal with respect to 
the Fermi-Dirac distribution at rest, 

w ^ = ePo/ktr + ! • ( 3 ) 

For the following derivations it is useful to choose nat- 
ural units, c = ks = h = 1. In addition, we will consider 
only the case for /i = 0, although a general approach is 
straightforward. By introducing a reference temperature 
T , we define 9 = T/T , p = p a /T , v = p/\p[, and using 
p° = \p\, we rewrite the equilibrium distribution as 

f eq , E {t,x,p,xf) = ePl(1 ^. s)/e + 1 , (4) 

where the subscript E stands for "Exact" . The distribu- 
tion f e q^E is expanded using tensorial polynomials P'"', 
for the angular contribution, and , for the radial de- 
pendence, such that 

1 00 

/ eg , E (i, = ^£af*>(t,f)pW(tf)F<*>(p) . 

7i, k 

(5) 

Here, the (2+l)-momentum vectors have been expressed 
in polar coordinates, p M = (p, p cos 0, p sin 4>) with v — 
(cos</>,sin</>) being a unit vector that carries the angular 
dependence </>, and the index i denotes a family of indices 
ii, i n € {1, 2} whose total number equals the order n 
of the tensor for the angular dependence, i.e. and 

af 1 ^ are tensors of rank n. Such an ansatz has been used 
by Romatschke et al. [TH] to expand the Maxwell- Jiittner 
distribution. Employing the Gram-Schmidt procedure, 
the radial polynomials F 1 -^ are constructed satisfying the 
orthogonality relation 

f™ ~™(P)F {k) (p)FW(p) = rP6 kl , (6) 

and the angular ones by satisfying 

£ T ^p^)pn {(j))=T {^p mn . (7) 

The resulting polynomials and T-constants up to second 
order are given in Appendix [X] With these polynomials 



3 



0.56 
0.54 
0.52 
■ 0.5 
0.48 
0.46 
0.44 







feq,E 
feq 



1 1 I 1 1 1 1 I 1 1 1 1 I 1 1 

0.5 1 1.5 2 



O/n 





O/n 



FIG. 1. Comparison betweem the expanded Fermi dirac 
distribution f eq and the full version feq.E as a function of the 
angular component <j>, for p = (left) and p — 3.5 (right), 
with = 1.0 (top) and 9 = 1.5 (bottom), u 1 = 0.0, u 2 = 0.05. 



and taking into account Eq. ^ , one can show that up to 
second order in n and k, we get 



(n/c) _ g 



(n) 



dp d0 

47T 27T 
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with = 1, #W = 2, and = 4. The explicit 
form of cS nk ^ is given in Appendix |b| Using Eq. the 
definitions of the polynomials, and their orthogonality 
relations it can be easily shown that the moments up to 



second order can be written in te 
aj with n, k < 2 (see Appendix 

truncated expansion of the distribution f eq up to second 
order becomes 



:ms of the coefficients 
BL and therefore, the 
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This is sufficient to recover the moments 
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= (e + P)U a U - Prf* 



(9) 
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of the full Fermi-Dirac distribution, Eq. Q. In Eq. ( 10b| , 
we have introduced the Minkowski metric tensor r/ ap , the 
particle density n = j$T 2 , the pressure P — ^p-nT and 
the energy density e = 2P, where £ denotes the Riemann 
zeta function, £(3) s» 1.202. 

Fig. [T] shows that the quality of the matching between 
the truncated f eq and the exact f e q,E, for P ~ 0, is very 



poor, in contrast with the case, p ~ 3.5. However, this is 
not surprising, since we are dealing with a gas of ultrarel- 
ativistic particles which are always moving at the Fermi 
speed, and therefore none of them has energy p = 0. On 
the other hand, the matching is reasonable for 6=1, 
while being off for 8 > 1. Thus, we conclude that 9 = 1 
offers the best approximation, and therefore, we will work 
with that value. In addition, we have found that 9 can- 
not be chosen far below unity because f eq can present 
negative values. The fact that 9=1 implies that the 
reference temperature T should be equal to the temper- 
ature of the electronic gas T. 



B. Momentum space discretization 

We now need to discretize the momentum space into a 
finite number N of discrete momentum vectors, p 1 ^ (with 
q = 0,...,N) such that we can replace integrals in the 
continuum momentum space by sums over a small num- 
ber of discrete momentum (2+l)-vectors. In order to do 
that, we use the Gaussian quadrature |29) . As an exam- 
ple, for the radial dependence of the expansion, in order 
to satisfy 
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(11) 

for k, I < 2, we should calculate the discrete p g > and 

(p) 

respective radial weights ui q , ■ By using the Gaussian 
quadrature theorem, we found the following values: 

pi = 0.484, u)f ] = 0.0369 
p 2 = 2.447, 4 ?) = 0-0176 
p 3 = 6.424, Jv } = 0.000719 . 

Note that in fact, p is always larger than zero, as ex- 
pected for ultrarelativistic particles, (see Appendix [C] for 
numerical values with higher precision). 

On the other hand, by following a similar procedure, 
we can calculate the N' discrete angles <fi q >> and angu- 
lar weights uj q t) (with q" = 1, ...,N'), such that, for the 
angular integrals over P( n '{vi) l {vj) m , one gets 
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P (B) K)'(»i) m = 



N' 
q"=0 



(12) 

where vi^u denotes Vi(q"). The above expression is rc- 



quired to be an exact quadrature formula for n < 2, 
and I + m < 2. The results for the discrete angles and 
weights functions are <j) q >> = f + (<?" — l)f and u q t) = g 
with N' = 6. 

By combining the radial and angular dependence 
of the discrete momentum (2+l)-vectors we get a to- 
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FIG. 2. The populations / q are moved between the nodes of 
a hexagonal lattice which are linked by the vector e^St. 



To(pq>,p q > COS ( 
the index q 



,Pq> 



where we have introduced 



(q',q"). This lattice cell configuration is 



shown in Fig. [2j where we can observe that for recov- 
ering hydrodynamics in graphene, we need a hexagonal 
lattice. This is a very convenient result, since due to the 
fact that it possesses the same honeycomb lattice sym- 
metries of graphene, we can reproduce with good accu- 
racy boundary conditions when modeling nanoribbons or 
other complex structures. 



The exact quadrature relations, Eqs. (11) and (12), 



ensure that the moments up to second order are still rep- 
resented exactly: 
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into a lattice, and each time step of 8t = 1 corresponds 
to one execution of the following steps: 

1. Calculate the equilibrium distributions f eq _q(t,x) 
from Eq. ^ using the macroscopic variables n — 
n(t,x), u = u(t,x), and T(t,x). At t = 0, n(t = 
0,x), T(t = 0,x), and u(t = 0,5?) are imposed as 
initial conditions. 

2. Collision: Introducing the post-collisional distribu- 
tions /q, calculate 



fq(t,x) = fq{t,x) 



p a U a 



n-iU&x) - f Eq ^(t,x)). 



At t = 0, take / q = / eg , q . 

3. Streaming: Move the / q along e q : 

/ q (t+l,f +e q ) = / q (t,i?) 

4. Calculate the new macroscopic variables. First we 
compute the energy density of the system by solv- 
ing the eigenvalue problem, (p a p l3 )U a = ell 13 , ac- 
cording to the Landau-Lifshitz decomposition |19) . 
From this, we get e and U a . Next, we use the rela- 
tion n = (p a }U a — n to obtain the particle density. 
Here, the average values, (p a ) and (p a p^}, are sim- 
ply 

q W (Pq) q 



q 
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(p°p 
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We have expanded and discretized the Fermi-Dirac 
equilibrium distribution for ultrarelativistic particles. 
Now, we will proceed to discretize the Boltzmann 
equation and find the evolution equation for the non- 
equilibrium distribution. 



C. Lattice Boltzmann algorithm 

With the expanded distribution functions and the dis- 
cretization of momentum space at hand, we may use the 
following discrete Boltzmann equation [TH [TH1 [52] , 



/ q (t+£t,x+e q $*)-/ q (t,a?) = - 
where we have introduced the notations e 



p a U a 

—p^riflfa &)-feq,<l(t, X)), 



(14) 

q = Pq/P°, 



and fq(t,x) = f(t,x,pq). Note that e q are unit vec- 
tors, which means that there are effectively 6 different 
e q . The discrete Boltzmann equation is now embedded 



The streaming step indicates that if we discretize the real 
space based on a hexagonal lattice where the sites are 
linked by e^St, as shown in Fig. [2j the values of / q will 
be moved between these sites exactly. This is known as 
"exact streaming" and crucial for the computational ef- 
ficiency and accuracy of the lattice Boltzmann methods, 
because it removes any spurious numerical diffusivity. 

In summary, we have developed a (2+l)-dimensional 
relativistic lattice Boltzmann scheme with the remark- 
able feature that it takes into account the Fermi-Dirac 
statistics, while recovering all the moments up to sec- 
ond order. The discretization is realized on a hexagonal 
lattice such that exact streaming is achieved. The fact 
that the quadrature corresponds to a hexagonal lattice 
allows to represent complex boundaries more precisely 
in graphene applications. This will be studied in more 
details in future works. 

Up to now, we are working with undoped graphene, 
/j = 0. However, by using the same orthogonal polyno- 
mials, we can easily integrate the Fermi-Dirac statistics 
for the doped case, obtaining the extended formulation. 
In this work, we will use [i = 0, in order to compare the 
results with previous models in the literature that use 
the Maxwell- Jiittner distribution, since transport theory 
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shows that in the case of undoped Fermi-Dirac statis- 
tics, the transport coefficients, namely shear viscosity 
and thermal conductivity, have the same expresions than 
for the Boltzmann statistics [TU]. Therefore, the shear 
viscosity takes the value of rj = (3/5)P(r — St/2) [3"U] . 
Later, we will use the doped case to study the Riemann 
problem, which to best of our knowledge has never been 
studied before. However, it is present when, for instance, 
laser beams are pointed to the graphene sheet in order 
to measure transport coefficients [3Tj . 



III. VALIDATION: RIEMANN PROBLEM 

In order to validate our model, we solve the Riemann 
problem for the ultrarelativistic Fermi-Dirac gas. The 
Riemann problem is a standard test for both, relativistic 
and non-relativistic hydrodynamics numerical schemes, 
because it involves the evolution of two states of the fluid 
initially separated by a discontinuity. In our case, we set 
up an effectively one-dimensional system of L x x L y = 
3000 x 2 nodes, using periodic boundary conditions in 
x and y components. Initially, there are two regions 
with particle densities, n = 1 {3L X /A > x > L x /4), 
and «i = 0.41 (x < L x /A and x > ZL X /A) creating a 
rectangular plateau of non-zero particle density in the 
center of the simulation zone. Here we consider and ini- 
tial constant temperature, To = 1. The initial velocity 
is set to zero and the value of the relaxation time r is 
calculated for two different values of £ = r]/(PoSt), with 

Pq = 9l * a 3 tioTq. The evolution of the system is displayed 
in Fig. [3] after 470 time steps, showing the generated 
shock wave. We have only plotted the region x > L x /2 
since the other one does not give additional information. 
Note that there is excellent agreement with the solutions 
provided by the model proposed in Ref. [24] for the same 
initial conditions. 
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FIG. 3. Density, pressure and velocity profile for the solution 
of the Riemann problem. Here £ = rj/(Po5t) is a dimension- 
less number. The expected results were calculated using the 
model in Ref. [24]. 



Let us now consider the case when the chemical poten- 
tial \x is different from zero. For this purpose, we follow 
the same procedure described before but this time, we 
keep /i 7^ 0. The development is straightforward, and 
therefore does not deserve a full explanation. The poly- 
nomials are the same as described in Appendix [XJ and 
the coefficients a[ nk ^ are calculated by using Eq. j|J. 

The hydrodynamic approach of electrons in graphene 
works for low doping, ^t/fceT <C 1 [THS]- Therefore, we 
can expand the discrete equilibrium distribution in pow- 
ers of n/ksT up to third order, neglecting errors of the 
order of (/i//csT) 4 . We perform additional simulations of 
the Riemann problem with the same parameters as be- 
fore, but now, varying the chemical potential. As we can 
observe from Fig. |2J increasing the chemical potential 
tends to increase also the effective viscosity of the sys- 
tem, smoothing the profiles of the velocity, pressure and 



density. This result is very interesting because it sug- 
gest that, in fact, impurities with soft potentials (small 
fi/ksT) in graphene samples can be treated as local mod- 
ifications in the effective viscosity of the electronic fluid. 
In other words, this result suggests a promissing way to 
include impurities in the hydrodynamic approach of elec- 
trons in graphene. Note that in this figure, there is a 
noise in the profile of the particle density. This numer- 
ical instability remains with the same amplitude and is 
always located at the boundary when n = no, and there- 
fore, it does not destroy the stability of the simulation. It 
can be due to the relevance of higher order terms which 
are not recovered by our expansion. 
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FIG. 4. Density, pressure and velocity profile of the solution 
of the Riemann problem, for different values of the chemical 
potential fi. 



are needed to recover the appropriate hydrodynamics. 
However, higher order moments can also be recovered by 
using the same procedure in this paper. 

In addition, we have developed a new lattice kinetic 
scheme to study the dynamics of the electronic flow in 
graphene. The model is validated on the Riemann prob- 
lem, which is one of the most challenging tests in nu- 
merical hydrodynamics, presenting excellent agreement 
with previous models in the literature. By increasing the 
chemical potential, we have found that the profiles of ve- 
locity, particle density, and pressure, change similar to 
the case when the viscosity is increased, concluding that 
increasing the Fermi energy results in increasing the effec- 
tive viscosity of the electronic fluid. This result suggests 
that soft impurities in graphene samples can be treated 
as local modifications of the viscosity, however, further 
studies must be performed in order to confirm this state- 
ment. 

The fact that we can propagate the information from 
one site to another in an exact way, avoiding interpola- 
tion, removes any kind of spurious numerical diffusivity. 
Therefore, we expect this model to be appropriated to 
study many problems in electronic transport in graphene 
in the framework of the hydrodynamic approach, e.g. 
turbulence and hydrodynamical instabilities in graphene 
flow, just to name a few. 

Extensions of the present model to take into account 
higher order moments of the Fermi-Dirac equilibrium dis- 
tribution as well as the inclusion of the distribution and 
dynamics of holes, will be a subject of future research. 
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Appendix A: Polynomials and T-constants 



V. CONCLUSIONS 

We have derived a new family of orthogonal polynomi- 
als using as weight function the Fermi-Dirac distribution 
for ultrarelativistic particles in two dimensions. By ap- 
plying the Gaussian quadrature we have calculated the 
set of representative momentum (2+l)-vectors, which al- 
lows us to replace the integrals over the continuum mo- 
mentum space by sums over such vectors. As a very 
interesting result, we have found that those vectors pos- 
sess the same symmetries than the honeycomb lattice 
of carbon atoms in graphene, making possible the accu- 
rate implementation of complex boundary conditions in 
future applications, such as point defects and nanorib- 
bons. The derivation has been performed by imposing 
that the expanded distribution should fulfill at least the 
first three moments of the equilibrium distribution, which 



In this section, we write explicitly the family of polyno- 
mials, which are orthogonal using as weighting function 
the Fermi-Dirac distribution at rest, with their respec- 
tive normalization factors. For the case of the angular 
dependence, we have 

P(°\v) = 1 

Pif(v) = ViVj - -Sij 
with normalization factors, 

= i 

r (1) =-s- 

1 p,ij - 2 13 
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For the case of the radial dependence, we have the 
polynomials 

F<°>(p) = 1 

fW(p)=p-c w , 



and, 



cio 



121og(2) ' 
F (2) (p) = p 2 - c 2 ip - c 20 , 

6(77r 4 log(2) - 15tt 2 C(3)) 



C21 



C20 



5(7r4-2161og(2)C(3)) ' 
7tt 6 - 3240C(3) 2 
10(7r 4 -2161og(2)C(3)) ' 
with C denoting the Riemann zeta function. The normal- 
ization factors for these polynomials are: 

(0) log(2) (1) ir 3 3C(3) 

F ~ in ' F 5761og(2) + 8^ ' 
r( 2) _ 1 fi9n s log(2) - 210tt 6 C(3) + 48600C(3) 3 ) 



400ttV tt 4 - 2161og(2)C(3) 

2250C(5) ' 



Appendix B: Coefficients for the expansion of f eq 
and relation to moments 



by 



The coefficients of the expansion in Eq. (|9| are given 
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7+1 

af ] = [( 2 " W(7 + 2) - cry] [ 7 4 (^, - ^/2) 

+ %/2] , 

a (° 2 ) = a( 2 )0 [/3 (2|1) (0 7 - 1) + /? (2|2) ((3^7 - 2) - 6) 
+ /3( 2 l 3 >(0 2 (3 7 2 -l)-2)] , 

aC«) = ^^) W7 + l)-l) 

+ /3( 2 l 2 ^(30 7 -2)( 7 + l) 

+ / 3( 2 l 3 )(30 2 7(7 + l)-2)]z il 7 , 

a i 22) = ^ (2|1) (( 2 - + 2 ) - ( 2<5 « + 

+ /3( 2 I 2 )(3# 2 ( 7 + l) 2 - 2(2 - %)(0( 7 + 2) - 
+ /3( 2 I 3 )(30 2 (7 + l) 2 - 2tr tf )] [ 7 2 (u^ - i^) 7 2 
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where cry = {—lyz.w.j or 



KO = ( 1 -1 ) ' 



(Bl) 



12tt 2 log(2) 



210tt 6 C(3) 



2161og(2)C(3) -tt 4 ' 
= 5 [2250C(5)(2161og(2)C(3) - tt 4 ) 
- 49tt 8 log(2) - 48600C(3) 3 ]^ 1 , 

0<?\ l ) = -147r 6 log(2) , 

/?( 2 I 2 ) = -15tt 4 C(3) 2 , 

^(2|3) = 3 2401og(2)C(3) 2 , 



which are approximately, a 1 - 1 ' w 0.994, ( 0( 2 I 1 ) w 
-1.629, a^ 2 ! 2 ) « -0.307, and a^ 2 ! 1 ) sa 0.567. 

To obtain the moments from the expansion of f eq , we 
expressed them in terms of the af 1 ^ using Eqs. ((§, 
([9]), and the expressions in Appendix [Aj e.g. (p°) — 
^(rWaW+c^ot")). 

Note that for the calculation of the coefficients a- 
we should use of the integration formula 



„n-l 



dx- 



-z- l a- n T{n)U n {-z) , 



which holds for n > 0, a G K, a > 0. Here, T(n) denotes 
the gamma function, which becomes T(n) = (n — 1)! for 
n£N. Li„(z) is the polylogarithm which can be defined 

using a power series: Li„(z) = J2T=i fk- If we consider 
the chemical potential in the Fermi-Dirac distribution to 
be zero, we have z = 1 and the relevant values of the 
polylogarithm become Lii(— 1) = — log(2), Li 2 (— 1) = 
-f^, Li 3 (-1) = -|C(3). On the other hand, for fj, ^ 0, 



we take z 



Appendix C: Results for radial Gaussian quadrature 

When the radial Gaussian quadrature is applied, the 
following values for the discrete pi are obtained: 

pi = 0.4840534751554060637550794361591 , 
p 2 = 2.4467448689670852668751189804200 , 
p 3 = 6.4243522612255152565859012563254 , 

with its respective weight functions 

lo[ p) = 0.0368730611359638360101542425978 , 
= 0.0175666801777458993453757617390 , 
J 3 P) = 0.0007191587244531629935841036927 . 



K. Novoselov, A. Geim, S. Morozov, D. Jiang, M. Kat- 
snelson, I. Grigorieva, and S. Dubonos, Nature Letters 
438, 197 (2005)| 

K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, 
Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. 
Firsov, |Science 306, 666 (2004)] 

D. P. DiVincenzo and E. J. Mele, Physical Review B 29, 
1685 (1984). 

S. Das Sarma, S. Adam, E. H. Hwang, and E. Rossi, 
|Rev. Mod. Phys. 83, 407 (2011)| 

A. K. Geim and A. H. MacDonald, Phys. Today , 35 
(2007). 

M. Miiller, J. Schmalian, and L. Fritz, Physical Review 
Letters 103, 025301 (2009). 

M. Miiller and S. Sachdev, |Phys. Rev. B 78, 115419| 
(2008) 



L. Fritz, J. Schmalian, M. Miiller, and S. Sachdev, Phys. 
Rev. B 78, 085416 (2008)" 



M. Miiller, L. Fritz, and S. Sachdev, [Phys. Rev. B 78 



115406 (2008) 

M. Miiller, J. Schmalian, and L. Fritz, Phys. Rev. Lett 



103, 025301 (2009) 



M. Mendoza, H. J. Herrmann, and S. Succi, Phys. Rev 
Lett. 106, 156601 (2011)] 



M. Mendoza, H. J. Herrmann, and S. Succi, [Sci. Rep. 3 
1052 (2013)1 

M. Mendoza, B. M. Boghosian, H. J. Herrmann, and 

S. Succi, Phys. Rev. Lett. 105, 014502 (2010). 

M. Mendoza, B. M. Boghosian, H. J. Herrmann, and 

S. Succi, |Phys. Rev. D 82, 105008 (2010)| 

D. Hupp, M. Mendoza, I. Bouras, S. Succi, and H. J. 

Herrmann, |Phys. Rev. D 84, 125015 (2011)| 



[16 
[17 



[18 

[19; 



[20 
[21 



[22 
[23 

[24 
[25 

£ 

[27 
[28 
[29 
[30 
[31 



Z. Xu and C. Greiner, Phys. Rev. C 71, 064901 (2005). 
I. Bouras, E. Molnar, H. Niemi, Z. Xu, A. El, O. Fochler, 
C. Greiner, and D. H. Rischke, Phys. Rev. Lett. 103, 
032301 (2009). 

P. Romatschke, M. Mendoza, and S. Succi, |Phys. Rev.| 
C 84, 034903 (2011)| 

C. Cercignani and G. M. Kremer, The Relativistic 
Boltzmann Equation: Theory and Applications (Boston; 

Basel; Berlin: Birkha user, 2002). 

X. He and L.-S. Luo, |Phys. Rev. E 56 6811 (1997)| 

N. S. Marty s, X. Shan, and H. Chen, |Phys. Rev. E 58,| 

6855 (1998)| 

X. He and L.-S. Luo, Physical Review E 55, R6333 
(1997). 

S. Succi, The lattice Boltzmann equation for fluid dynam- 
ics and beyond (Clarendon Press and Oxford University 
Press, Oxford and New York, 2001). 
M. Mendoza, I. Karlin, S. Succi, and H. J. Herrmann, 
|Phys. Rev. D 87, 065027 (2013)| 
M. Y. Han, B. Ozyilmaz, Y. Zhang, and P. 
206805 (2007)1 

Hod, and G. E. Scuseria, Nano Letters 6, 



Kim, 



Rev. Let t. 98, 
V. Barone, O. 
2748 (2006) | 

Z 



Phys. 



Physik (Zeitschrift fur Physik) 47, 542 



F. Jiittner, 
(1928). 

H. Grad, Communications on Pure and Applied Mathe- 
matics 2, 331 (1949). 

P. J. Davis and P. Rabinowitz, Methods of numerical in- 
tegration, 2nd ed. (Academic Press, Orlando, 1984). 
M. Mendoza, I. Karlin, S. Succi, and H. J. Herrmann, 
J STAT 2013, P02036 (2013). 

J.-U. Lee, D. Yoon, H. Kim, S. W. Lee, and H. Cheong, 
Phys. Rev. B 83, 081419 (2011)| 



